In-class Exercise 5

This in-class exercise is a revision of what is covered in Hands-on Exercise 4, which revolves around spatial point patterns analysis.

Installing and Loading Packages

• maptools: creating spatial* objects
• sf: importing geospatial data
• raster: creating raster object
• spatstat: performing SPPA
• tmap: plotting thematic maps
• tidyverse: data handling
• plotly: creating interactive plots
• ggthemes: provides more advanced themes for ggplot

packages = c('sf', 'spatstat', 'raster', 'maptools', 'tmap', 'tidyverse', 'plotly', 'ggthemes')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}

Importing Geospatial Data

• st_read() of **sf* package is for importing geospatial data; not appropriate for importing aspatial data.
• The output object is in tibble sf object class.

mpsz_sf <- st_read(dsn="data/shapefile",
                   layer="MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\Users\user\Desktop\SMU\Y4S1\IS415\geniceee\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

• Projection is in SVY21.

Importing Aspatial Data

read_rds() of readr package is used instead of readRDS() of base R is used. This is because the output of read_rds() is in tibble object.

childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")

Note: There are some data issues in childcare dataframe because Lat and Lng should be in numeric data type. The coordinate fields seem to be in decimal degrees and it is assumed that WGS84 referencing system is used.

Convert Aspatial Data into sf objects

CHAS_sf <- st_as_sf(CHAS, 
                    coords=c("X_COORDINATE", "Y_COORDINATE"),
                    crs=3414)

Note: st_as_sf() accepts coordinates in character data type

childcare_sf <- st_as_sf(childcare, 
                    coords=c("Lng", "Lat"),
                    crs=4326) %>%
  st_transform(crs=3414)

• No significant difference if the columns are not numeric during conversion. But we should convert the columns to numeric type for best practice.

Plotting for Reviewing

• alpha: set transparency of dots
• color: set colour of dots
• size: set size of dots

tmap_mode('view')
tm_shape(childcare_sf) + 
  tm_dots(alpha=0.4,
          col="blue",
          size= 0.05) +
tm_shape(CHAS_sf) + 
  tm_dots(alpha=0.4,
          col="red",
          size= 0.05)

Geospatial Data Wrangling

Converting from sf to Spatial* dataframe

Using as_Spatial() of sf package

• coordinates: geometries object
• bbox: extent of geometries
• proj4string: defines projection

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

Converting Spatial* dataframe into Spatial* objects

Using as.SpatialPoint() / as.SpatialPolygon() of maptools package

• Input (e.g. childcare) needs to be a Spatial* dataframe
• Output object is not a dataframe anymore since it does not contain data object

childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")

Converting Spatial* objects into ppp objects

Using as.ppp() function of maptools package

• Drops away projection information

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

Remove Duplicate Points using Jitter

childcare_ppp_jit <- rjitter(childcare_ppp,
                             retry=TRUE,
                             nsim=1,
                             drop=TRUE)

any(duplicated(childcare_ppp_jit))

[1] FALSE

CHAS_ppp_jit <- rjitter(CHAS_ppp,
                             retry=TRUE,
                             nsim=1,
                             drop=TRUE)

any(duplicated(CHAS_ppp_jit))

[1] FALSE

• Cannot use tmap() to plot ppp objects for visualisation. Need to convert the ppp object to sf object before it is possible to plot using tmap(). Will have to resort to using plot()

Extract Punggol Planning Area

• @: to access data object
• $: to access columns in data object

pg <- mpsz[mpsz@data$PLN_AREA_N=="PUNGGOL",]

Convert SpatialPolygonsDataFrame into Spatial Polygons object

pg_sp <- as(pg, "SpatialPolygons")

Convert SpatialPolygons into owin object

pg_owin <- as(pg_sp, "owin")

Extract Spatial Points within Owin

• Points will now within the Punggol sub-district

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]

plot(childcare_pg)

L-function

L_childcare <- envelope(childcare_pg,
                        Lest,
                        nsim=99,
                        rank=1,
                        global=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

L_CHAS <- envelope(CHAS_pg,
                        Lest,
                        nsim=99,
                        rank=1,
                        global=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_childcare)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
}

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_CHAS)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
}